Bloch-mode analysis for retrieving effective parameters of metamaterials 
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We introduce a new approach for retrieving effective parameters of metamaterials based on the 
Bloch-mode analysis of quasi-periodic composite structures. We demonstrate that, in the case of 
single-mode propagation, a complex effective refractive index can be assigned to the structure, being 
restored by our method with a high accuracy. We employ both surface and volume averaging of 
the electromagnetic fields of the dominating (fundamental) Bloch modes to determine the Bloch 
and wave impedances, respectively. We discuss how this method works for several characteristic 
examples, and demonstrate that this approach can be useful for retrieval of both material and wave 
effective parameters of a broad range of metamaterials. 
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I. INTRODUCTION 

The study of artificially structured metamaterials 
(MMs) attracts attention of scientists and engineers due 
to their unprecedented electromagnetic properties. Neg- 
ative refractive index, very large or near zero values of 
both permittivity and permeability, giant optical activ- 
ity - these are just a few examples of the properties which 
MMs can provide^. As was established, it is convenient 
to describe the MM properties by employing the concept 
of effective parameters (EPs), such as refractive index 
n, impedance z, permittivity e and permeability /x, pro- 
vided that these EPs can be introduced 2 -. The EPs sim- 
plify significantly the description of the MM properties, 
including the propagation of electromagnetic waves in- 
side a MM slab and their reflection and transmission at 
the MM boundaries. 

The state-of-the-art of homogenization infers that re- 
trieved EPs are of two types^r— : 

(i) Material (or local) effective parameters (MEP) em 
and hm- They give the relation of the field vectors 
D = £j\/£()E and B = ^mA'oH. The material effective 
parameters show the evolution of the wave inside a meta- 
material. Material EPs depend only on the properties of 
the material (we do not consider here the problem of the 
Drude transition layers^). Specifically, material EPs are 
important, for example, for the superlens performance of 
the slab with negative refractive index^. The relations to 
the refractive index n and wave impedance Zw are: 



z w 



(1) 
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(ii) Wave (or nonlocal) effective parameters (WEP) Ew 
and fj,w- They are usually restored from the reflection 
and transmission coefficients of a MM slab 6 being as- 
signed as the parameters of the corresponding homoge- 
nous slab. Sometimes this approach leads to violation 



of locality conditions, and this situation was actively dis- 
cussed in the literature 2 - - — . The WEP may allow one 
solving the scattering problem (reflection/transmission 
determination) for a MM slab of another thickness. They 
often depend on the thickness of the MM slab (in terms 
of the number of unit cells, see e.g. Ref. [Hj]), with only 
rare exceptions^. 

For a homogeneous medium with the structural ele- 
ment characteristic size a, which is much less than the 
wavelength A, the material and wave EPs are the same. 
However, in many practical cases MM's unit cell is only 
a ~ A/10— A/4 and material and wave parameters are not 
equivalent to each other—. It is obvious that the reflection 
from a MM slab should depend on whether the MM slab 
termination coincides with the border or with another 
cross-section somewhere in the middle of the unit cell, so 
the wave EPs depend on the MM opening cross-section. 

The knowledge of the WEP and MEP is needed for 
development of metamaterial based devices. This would 
be desirable to obtain both sets of EPs within a simi- 
lar simple calculation procedure. The importance of the 
EPs restoration is emphasized by a variety of the existing 
retrieval methods, which are summarized in the Table HI 



This paper aims to introduce and discuss in detail 
an approach described in [5, 38] for retrieving the wave 
and material effective parameters. First, we calculate 
the dispersion bands of the long enough periodic me- 
dia by employing the high-resolution spectral analysis 
method22r— . This method is developed for periodic 
structures (in fact, quasi-periodic, taking into account 
a finite size of the structure) composed of arbitrary unit 
cells. After defining the dispersion of the dominating 
Bloch modes, we introduce a complex refractive index 
which can be attributed to the effective parameters of 
the metamaterial with a high accuracy. 

Next, we introduce an effective impedance. Follow- 
ing Refs. 2.38], we apply the volume or surface averag- 
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TABLE I: Comparison of the EPs restoration methods 



Method and References. 


Effective parameters type 
and comments 


Reflection-transmission 
(Nicholson- Ross- 
Weir (NRW) )*&!*>!* 


WEP: Scalar, restored for 
normal or inclined incidence. 


Wave propagation 16 ' 17 


WEP: Scalar, restored for 
normal incidence. 


Field averaging**"-* 2 - 


MEP: Scalar or tensor. 


Analytical and semi- 

l i4, 23-26 

analytical*** — 


MEP: Tensor. 


Single interface scattering* 7 - 


WEP: Scalar, restored for 
normal incidence. 


Non-local dielectric 
function 2 ^ 


MEP: Nonlocal dielectric 
function, tensor. 


Current-driven^** 6 - 


MEP: Tensor. 


Quasi-mode* 7 - 


WEP: Scalar, restored for 
normal incidence. 



ing of the electric and magnetic fields of the dominating 
Bloch mode, which leads to the wave or Bloch (input) 
impedance EPs retrieval. Having both refractive index 
and impedance, we restore effective permittivity and per- 
meability accordingly to Eqs. ([TJ and ([2]), which will be 
either MEP or WEP, respectively. However, in contrast 
to the refractive index retrieving, the wave impedance 
retrieving procedure may encounter problems especially 
in application to MMs with the negative refractive in- 
dex. Caution should be paid to the fields computed 
via direct numerical solution of Maxwell's equation by 
Maxwell's solvers. For example, in the CST Microwave 
Studio, which we used, the returned magnetic field cal- 
culated on a grid is magnetic induction b//x and not 
magnetic strength h. Ignoring this fact when restoring 
the impedance from the electric and magnetic fields ratio 
can cause the real part of impedance becoming negative 
in the region of the negative refractive index, and cor- 
respondingly the negative energy flux is obtained. Such 
flux behavior is connected with its definition through the 
H field, the fact that was emphasized by Silverinha et 

.,128,32 



The paper is organized as follows. In Sec. [H] we for- 
mulate the general concept and technical details of our 
approach. The successful MEP retrieving examples in 
the case of homogeneous media and different types of 
composite MMs are summarized in Sec. IIIII In Sec. IIII1 
we also present the examples when the wave impedance 
retrieval leads to incorrect interpretation of EPs and, as 
a consequence, it connects impedance with the energy 
flux with wrong flux direction. Finally, in the concluding 
Sec. IIVI we discuss both advantages and constraints of 
the novel approach introduced here. 



II. GENERAL APPROACH 

The dispersion analysis is based on the Bloch modes 
expansion of the field propagating inside a MM slab. We 
simulate the field propagation by the commercial CST 
Microwave Studio software** with the finite-integrals 
Maxwell solver. 

We excite the MM slab, which consists of the period- 
ically arranged unit cells of the period a = (a x ,a yi a z ), 
with a plane wave propagating along the z— axis and elec- 
tric field polarized along x— axis (see Fig. [T]). In princi- 
ple, the slab may be arbitrarily thick, but not less than 
3-4 MM monolayers for that we can neglect the so-called 
Drudc transition layers 2 .. 

We use perfect electric, perfect magnetic and open 
boundary conditions for the x— , y— and z— boundaries 
respectively and the time-domain solver in calculations. 
A broadband Gaussian pulse is used as a field source. 
Only one simulation is needed for the whole spectrum 
calculation. The fields on different frequencies are calcu- 
lated through the Fourier transformations from the time- 
dependent signals collected with 3D field monitors. 

Let us consider the plane wave normally incident 
from vacuum onto the MM slab. Its electric E v = 
E v q exp(ifcoz) and magnetic H v — H v oexp(ikQz) fields 
are connected via the impedance of the free space, Zq = 
E v q/H v q — y/ /io/eo ~ 1207T Ohm. Here fco = oj/c- is 
the wavenumber of the free space, and we assume the 
exp(— icut) time dependence. 

In general case, several Bloch modes** - — may be ex- 
cited in the slab for each frequency oj, so the overall field 
may be represented as a sum 



A I 



E(r) = J2 E m(r), 



M 



H(r) = ]T H m (r), 



(3) 



(4) 



where m is the Bloch mode number, M is the total num- 
ber of excited modes, and r = (x, y, z). In the desirable 
case of local quasi-homogeneous MM there are only two 
modes in the slab: one forward and one backward propa- 
gating. A larger number of modes may be excited in the 
case of MM with strong spatial dispersion 2 . 

The field profiles of Bloch modes can be represented 



as***-** 



E m (r) = 



H m {v) 



p#0 



iGpz 



(•5) 



p#0 



„iK m z 



, (6) 
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FIG. 1: (Color online). Simulation configuration. Wave is 
normally incident from vacuum. Wave propagation and meta- 
material stacking direction is along z-axis. Electric field of the 
plane wave is polarized along x-axis. 



where K m is the Bloch wavenumber, G — 2ir/a z , p is an 
integer number. We note that the field representation in 
Eqs. ([3]), (|4]) is invariant with respect to a transformation 



K m — > K m + Gp' and E m>p — > E mtP+p > for an arbitrary 
integer p' . Accordingly, we can always select the value 
of K m such that E m fi is the largest harmonic amplitude, 
and we use this convention in the following. 

The key feature of the high-resolution spectral anal- 
ysis metho d 39 i 41 is decomposition of the total field ob- 
tained in simulations into a sum of Bloch modes, effec- 
tively inverting Eqs. <J3J) ^ (SI) - The only prior information 
required for the application of this method is the number 
of strongest Bloch modes excited in the structure (M). 
Then, through specialized numerical fitting described in 
Refs. 39.41] we extract wavenumbers K m and field pro- 
files E m (r) , H m (r) of all forward and backward propagat- 
ing Bloch modes at each frequency u). By monitoring the 
accuracy of such decomposition in terms of field match- 
ing, we check whether other ignored Bloch modes have 
significant excitation amplitudes, and if this is a case we 
increase the number M to take more modes into account 
and repeat the whole decomposition procedure. 

It is an important advantage of our approach that the 
standing wave, which is usually formed inside the slab 
due to the multiple reflections from the boundaries and 
brings the restrictions to the conventional wave propa- 
gation retrieval method^, is not an issue in the present 
case, since we can separate forward and backward propa- 
gating Bloch modes. In the following, we denote the field 
profiles of the dominant forward and backward waves as 

{E,H} + = {E,H} m+ , {E,H}. = {E,H} m _, (7) 

and the corresponding wavenumbers 



K, = K, 



K _ = K n 



(8) 



where m + and rn_ are the numbers of the dominant for- 
ward and backward Bloch modes, respectively. 

If several Bloch modes are excited and propagate in 
a MM, such composite cannot be homogenized and no 
meaningful EPs can be introduced. The homogeneity of 



MM and the influence of the higher-order Bloch modes 
have been discussed extensively in the Refs. [13U47H48| . 
However, if only one forward mode can be distinguished 
by the lowest damping, we can count it as the dominating 
one and neglect the presence of the higher-order modes. 
As a rule it is the fundamental Bloch mode. The numeri- 
cal criterion of homogeneity from the Bloch modes point 
of view was formulated in Ref. [Hj]. Another possibility 
to check the single mode regime is to calculate the mis- 
match 5 of the restored sum of forward and backward 
propagating fundamental mode fields, Ef = E + + 
and the original field E taken directly from numerical 
simulations: 



S = 



f\E- E f \ 2 dxdydz 
J \E\ 2 dxdydz ' 



(9) 



where integration is performed over the computation do- 
main. In all the case studies presented below the mis- 
match S is below 1.5%. So, in this manuscript we con- 
sider the MMs that have a dominant fundamental mode, 
and the higher-order Bloch modes can be neglected. Ac- 
cording to the concept of homogenization, we aim to 
find effective parameters for an equivalent homogeneous 
medium, where the wave propagation would be essen- 
tially the same as in the periodic structure. After deter- 
mining the propagation constant K + of the fundamental 
mode we assign our structured material with the effective 
refractive index n = K + /ko. 

The second part in restoration is connected with the 
effective impedance. We use the fields E + , H + of the fun- 
damental Bloch mode in the both Bloch zb and wave z\v 
impedances restoration. First, we perform fields surface 
averaging at the (x, y) cross-section of the simulated slab: 



Esa{z) = / E + (x,y,z)dxdy/a x a y , 
Js 



Hsa(z) = / H + (x,y, 
Js 



(10) 



z)dxdy /a x a y . (11) 



Taking the values of the fields E$A,j — EsA(zj), Hsaj — 
HsA(zj) at the unit cell borders Zj — ja z , where j is an 
integer number, we determine the Bloch impedance^: 



zb 



E SA,j 

ZqHsa. 



(12) 



Note that Bloch impedance zb does not depend on j, 
which can be checked by substituting Eqs. ([5]) and © 
into Eqs. {TDJ and (fl~Tj) . 

In order to restore wave impedance (zw), we need to 
calculate the volume averaged fields 2 i?vA and H\ 



zva, 



zw 



Eva 
ZqHva 



(13) 



Since the wavenumbers in the periodically structured 
and equivalent homogeneous media are equal, we need 
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to establish the correspondence of the field amplitudes 
in front of the common exp(iK + z) multiplier. Accord- 
ingly, we define the volume-averaged fields by perform- 
ing integration over a single unit cell with the multiplier 
exp(— iK + z) to cancel the phase evolution: 



E VA = E SA (z)exp(-iK + z)dz/a z , 

J z b 



(14) 



VA 



H SA (z)exp(-iK + z)dz/a z , (15) 



where z& is an arbitrary location inside the structure. We 
can also express the averaged fields through the harmonic 
amplitudes by substituting Eqs. ([5]) and ([6]) into Eqs. (fT4|) 
and ([TBI . 

E VA = / E m+fi (x,y)dxdy/a x a y , (16) 



H YA = / H m+:0 (x,y)dxdy/a x a y . (17) 



We see that the volume-averaged fields do not depend 
on Zb, as their values are defined through the dominant 
Bloch-wave harmonic amplitude which is z-independent. 

For the extraction of E and H fields from the CST Mi- 
crowave Studio simulations we use electric and magnetic 
field monitors. However, the raw microscopic magnetic 
field that CST returns is not h(r), but rather b(r)//xo as 
a straightforward solution of microscopic Maxwell's equa- 
tions. As this was shown by M. Silveirinha et al . 28 ' 32 , the 
employment of the volume averaged magnetic induction 
By^(r) instead of Hy^(r) can give an incorrect direction 
of the Poynting vector for negative index metamaterials. 

For the correct determination of the volume averaged 
magnetic field we employ the definition 



H 



B 



VA 



VA 



/'() 



-M 



VA, 



M 



VA 



Ut x 3)dV 



2V 



(18) 



(19) 



where M^^ - is the volume averaged magnetization vec- 
tor and J is the current density. In principle the mag- 
netization can be calculated by a numerical integration 
routine directly from the definition. However, we choose 
another, more elegant, approach following the findings of 
Silveirinha for the transverse-averaged magnetic fields 3 ^. 
First, we decompose My^ into two parts: along the di- 
rection of propagation (unit vector u z ) and orthogonal 
to it 



H 



B 



VA 



VA 



Mo 



- {M VA • u z )u z + u 2 x {u z x M VA ). (20) 



Then, we project the previous expression onto the tan- 
gential plane. Taking into account that the magnetic 



field has dominating polarization in the tangential plane 
provides 



H 



B 



VA 



VA 



Mo 



u z x (u z x M VA ) 



B 



SA 



Mo 



(21) 



This equation holds for the long-wavelength limi t 28 ' 32 . 
Thus, in order to calculate the correct values of the wave 
impedance (and Poynting vector) one can use volume av- 
eraged numerical electric field Eva, but surface averaged 
numerical magnetic field B$a 



z w = 



Eva^o 
ZqBsa 



(22) 



We would like to remark that Eq. ([22|) makes a bridge 
between our approach and that of papers with averaging 
field procedures^—, where effective magnetic functions 
are obtained via volume averaging of B fields, but surface 
averaging of H fields. 

Deriving effective permittivity and permeability from 
Eqs. Q, © we find the MEP of the metamaterial. Ac- 
cordingly reversing Eqs. ([T]), (|2|) for the Bloch impedance 
(fT5)) we ends with the set of metamaterial WEP. Thus, 



and 



£m = n/zw, t*M = nz w , 



e w = n/z B , fi w = nz B . 



(23) 



(24) 



The latter should be equal to these given by the NRW 
method^. We emphasize that determination of the propa- 
gation constants and impedances is straightforward, does 
not involve any inverse functions and is made on the ba- 
sis of the same simulated fields for both wave and Bloch 
impedances. 

We should mention a practical issue important for 
the implementation of the proposed approach. Comput- 
ing fields by the finite-difference or finite-integral time- 
domain methods we should take into account a phase 
shift between the electrical and magnetic fields connected 
with the staggered Yee mesh. The electric and magnetic 
fields are calculated at different time moments shifted by 
At/2, where At is the simulation time step. For the case 
of CST Microwave Studio, which we are using, the mag- 
netic field phase is always shifted by Ac/) = usAt/2, so 
we corrected the magnetic field values by corresponding 
phase factor exp(zo;At/2). 



III. SPECIFIC EXAMPLES OF 
METAMATERIAL STRUCTURES 

We tested our approach on several examples, start- 
ing with the simplest ones. The unit cells sketches of 
the designs are shown in Fig. [2] We considered: (1) ho- 
mogeneous slab [see Fig. HJa)] - two cases: lossless and 
Lorentz dispersion in s and /j, with negative index of re- 
fraction, (2) a set of the nanospheres with the plasmonic 
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FIG. 2: (Color online). Sketches of the materials designs con- 
sidered: homogeneous material (a), plasmonic nanospheres 
(b), split cube MM (c), wire medium (d), fishnet MM (e) and 
split cube in carcass MM (f). 

resonances [see Fig. Hfb)], (3) split cubes MM that pos- 
sess magnetic resonance and negative permeability [see 
Fig. [He)], (4) wire medium that gives negative permit- 
tivity [see Fig.[2ld)], (5) negative refractive index fishnet 
MM [see Fig.[3fe)], and (6) split cube in carcass MM [see 
Fig. [210] ■ I R all cases, the MM slab consisted of 10 mono- 
layers. For comparison, WEP for three-monolayers-thick 
slabs were calculated with the NRW method 6 . 



A. Homogeneous materials 

A slab of homogeneous material is the simplest object 
to test the retrieval approach, since the restored EPs can 
be compared with the exact values. 

A homogeneous slab was artificially divided into 10 
meta-atoms of the size a x — a y — a z = 100 /im. For the 
case of the homogeneous medium, the material and wave 



parameters are identical, so we should only compare the 
given constitutive parameters with the retrieved MEP. 

For the homogeneous lossless slab with constant pa- 
rameters: e — 4 and fj, = 1 the EPs were in a perfect 
agreement with the theoretical permittivity and perme- 
ability (not shown). The relative retrieval error was less 
than 0.2%, which can be attributed to numerical disper- 
sion effect in finite-difference numerical simulations. 

In another example, we consider the frequency dis- 
persive permittivity and permeability described by the 
Lorentz model: 

e(w) = £co + Estat—, — J ' ( 25 ) 

W 0e _ *7e^ - & 



/iM = /ioo + ton — o> ( 26 ) 

where =1, e sta t =1-7, uj 0e — 2tt x 198 x 10 9 s -1 , 

7e = 2-7T X 10 10 S _1 , ^oo =1, Vstat =1.3, UJ 0m = 2n X 

202 x 10 9 s" 1 , 7 m = 2tt x 10 10 s _1 . 

The restored effective parameters are in good cor- 
respondence with the original EPs [see Fig. [3]. The 
small differences are observed only in the resonant region 
around 200 THz where losses are high. The retrieval re- 
sults in Fig. [3] show that retrieving through the Bloch 
mode analysis is applicable to a range of materials with 
or without losses with positive and negative n, e and fi. 



B. Metamaterial composed of plasmonic 
nanospheres 

Metallic nanospheres possess plasmonic resonances. 
Being arranged in the regular structure, the nanospheres 
with a radius r«A make a MM. It is expected that the 
nanospheres MM should have the permittivity which is 
different from the host permittivity and its permeabil- 
ity should be close to 1, since the nanospheres are non- 
magnetic. 

The silver nanospheres of the radius r — 30 nm were 
placed in vacuum in the cubic array with the period a x = 
o>y = a>z — 200 nm. Silver was considered as the Drude 
metalSi with the plasma frequency u> p — 1.37 x 10 16 s _1 
and collision frequency 7 C = 8.5 x 10 13 s _1 (see Ref. [Hoj ) . 
The sketch of the design is shown in Fig. [5Jb). 

Effective refractive indices restored with the NRW 
method and our approaches are identical [see Fig.@|a)] as 
it was expected. Bloch impedance zb, retrieved with the 
field surface averaging [see Fig. Hp)), triangles] is identi- 
cal to the one restored with the NRW method, [Fig.QJb), 
solid lines]. 

There is a little difference between wave impedance zw 
[see Fig. H£b), circles] and zb (triangles). They experi- 
ence slight oscillations around the value of zw — 1 + Oi. 
As a consequence of that both permittivities exhibit res- 
onances around 660 THz, 690 THz and 730 THz [see 
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FIG. 3: (Color online). Retrieved effective parameters 
(circles) of the homogeneous medium with Lorentz disper- 
sion in permittivity and permeability: refractive index (a), 
impedance (b), permittivity (c) and permeability (d), real 
(black) and imaginary (green/grey) parts. Results are com- 
pared with the original values (solid lines). 

Fig. Etc)], but of different strength. At the same fre- 
quencies, the magnetic permeability shows non-physical 
negative imaginary part, so-called antiresonance behav- 
ior that normally would correspond to the gain in the 
system. However, material EPs Sm and fiM, restored 
via the volume-averaged fields are free from the antires- 
onances on frequencies up to 700 THz. Small negative 
values of 9(£m) are due to the calculation errors with 
the staircase approximation of the spherical shapes. 

The permeability which is supposed to be around 

1 since the nanospheres are non-magnetic, is indeed 
around 1 on frequencies up to 700 THz, but starts 
to oscillate on higher frequencies, especially at around 
750 THz [see Fig. Hfd)]. It looks as we have strong mag- 
netism from the non-magnetic MM consisting of electric 
dipoles. In fact, at frequency 750 THz the condition for 
the first Bragg resonance is satisfied, so the MM cannot 
be considered as homogeneous and cannot be assigned 
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FIG. 4: (Color online). Effective parameters of the MM 
consisting of plasmonic nanospheres: refractive index (a), 
impedance (b), permittivity (c) and permeability (d), real 
(black) and imaginary (green/grey) parts. Retrieved results 
by volume-averaged (circles) and surface- averaged (triangles) 
fields are compared with the NRW method (solid lines). 

with meaningful effective parameters — . 



C. Split-cube metamaterial 

We choose a split cube MM as an example of a mag- 
netic material with negative permeability in the infrared 
rang o 13 i 51 . The sketch of the design, which is a 3D 
generalization of the symmetric split-ring resonator—, is 
shown in Fig. [2jc). The cubic unit cell of a x = a y = 
a z = 250 nm consists of the silver thin-wall structures 
(Drude metal) embedded in silica (permittivity 2.25). 
The geometrical parameters were taken the same as in 
the Ref. [3. 

In the line with the previous cases, the refractive in- 
dices retrieved with different methods coincide, showing 
a resonance around 160 THz [see Fig. EJa)]. Bloch and 
wave impedances exhibit strong resonance behavior in 
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FIG. 5: (Color online). Effective parameters of the split cube 
magnetic MM: refractive index (a), impedance (b), permit- 
tivity (c) and permeability (d), real (black) and imaginary 
(green/grey) parts. Results by volume-averaged (circles) and 
surface averaged (triangles) approaches are compared with 
the NRW method (solid lines). 



FIG. 6: (Color online). Effective parameters of the wire 
medium: refractive index (a), impedance (b), permittivity (c) 
and permeability (d), real (black) and imaginary (green/grey) 
parts. Results by the volume- averaged (circles) and surface 
averaged (triangles) approaches and NRW method (solid line) 
are compared with the analytical predictions (stars). 



the area around 160 THz. A small peak in the impedance 
restored with the NRW method only at the frequency 
91 THz appears at the Fabry-Perot resonance of the slab 
and is a numerical artifact intrinsic to the S-parameter 
method [see Fig. EJb)]. The spurious peaks in the EPs 
due to Fabry-Perot resonances can be avoided with wave 
propagation methods as it was reported in Ref. [13] • 

Effective parameters restored via surface and volume- 
averaged fields expose strong antiresonance behavior for 
the effective dielectric permittivity. Such behavior ordi- 
nary for WEP cannot be accepted in assigned MEP. The 
reasons for very similar appearance of effective parame- 
ters revealed by formulas (|2"T) . (fT3|) we assign to a strong 
magnetic resonance, which brings domination of mag- 
netic field performance through B$a denominator and 
thus to formal equivalence of effective impedances. How- 
ever, the full picture of failure of formula (fT5)) has yet to 
be understood. 



D. Wire-medium structure 

Wire medium 53 is a well-known example of the 
negative-permittivity MM. In the case of the square lat- 
tice of perfectly conducting wires in vacuum, when radius 
of the wires r is much less than the unit cell size, r <C a, 
an analytical expression for the effective permittivity is 
given in Refj 5 ^: 

, . 2ttc 2 
£e//H = 1 " a ^(log^ + 0.5275)- (27) 

We simulated the r — 5 /im-radius wires made from 
the perfect electric conductor arranged in a square lattice 
with a x = a y = 500 fim in vacuum [see the sketch in 
Fig. HJd)]. Comparison of the retrieved and analytical 
EPs is presented in Fig. [5] 

Due to the rectangular spatial discretization of the 
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round-shaped wires in the simulations we see the differ- 
ence in the effective impedances retrieved through the 
field averaging procedure (both of them!) and the NRW 
method. It causes deviations in effective permittivities. 
Permittivity retrieved with the NRW method is closer 
to the analytical results [see Fig. Etc)]. What concerns 
permeability, the NRW method retrieves paramagnetic 
$t([iw) ~ 1-2 [see Fig. |5[d)], while the wire medium is 
expected to be a non-magnetic MM. Within the field- 
averaging approach the retrieved jiw perfectly coincides 
with the theoretical prediction, while hm seems to be 
more sensitive for the staircase approximation errors. 

We should note that because we study wave propaga- 
tion perpendicular to the wires no any spatial dispersion 
effect showed up during the restoration, and results are 
physically sensible. 



E. Fishnet metamaterial 

The fishnet MM^£ is one of the most promising 
negative-index metamaterials for the optical and infrared 
regions. It consists of the metallic double wires extending 
in the x— and y— directions [see the sketch in Fig.^e)]. 

We use the geometrical and material parameters of the 
fishnet MM from Ref. except adjusting the period in 
z— direction to a z = 150 nm. The unit cell transverse 
sizes are a x = a y — 600 nm. Silver layers (silver treated 
as the Drude metal) of the thickness 45 nm are separated 
with the MgFi dielectric with refractive index n = 1.38 
and thickness 30 nm. This metal-dielectric-metal sand- 
wich is placed in vacuum. 

The refractive indices retrieved with our approach and 
the NRW method are slightly different [see Fig. |7fa)]. 
This is not surprising since the NRW method is applied 
to a three monolayers-thick slab. It is well known that 
the thin-slab effective refractive index of the fishnet con- 
verges slowly to the bulk values with the increase of the 
slab thicknes o 12 ' 55 . Our approach based on field propa- 
gation in 10 layers gives the refractive index close to its 
bulk values. Bloch and wave impedances are different as 
well [see Fig.^b)]. We also expect that the NRM results 
would converge to ours if ten layers will be considered. 

Effective parameters obtained by both types of field 
averaging are quite close to each other. The feature of 
the fishnet behavior is the negative index region free from 
the anti-resonances both in e and /i. The NRW results 
exhibit hardly visible anti- resonance for 9(e), which is 
corrected via field-averaging procedures. 



F. Split-cube-in-carcass metamaterial 

A fishnet MM is an example of a medium with a neg- 
ative refractive index. To check that we can assign ef- 
fective parameters, which will not show any antireso- 
nances, we consider another negative-index metamate- 
rial with strong spatial dispersion, namely split cube in 
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FIG. 7: (Color online). Effective parameters of the fishnet 
negative- index MM: refractive index (a), impedance (b), per- 
mittivity (c) and permeability (d), real (black) and imaginary 
(green/grey) parts. Results by volume-averaged (circles) and 
surface averaged (triangles) approaches are compared with 
the NRW method (solid lines). 



carcas g 13 i 17 [see the sketch in Fig. HIT)]. Its remarkable 
property is extreme fast convergence of parameters such 
that its effective refractive index is the same for the 1- 
layer thick slab and for the bulk MM represented by 
the infinite number of layers. However, as was shown 
in Ref. [48j |. even being 3D cubic symmetric by design, 
split cube in carcass is anisotropic in the resonant region. 

The cubic unit cell of a x — a y — a z = 250 nm 
[Fig- Elf)] consists of the silver split cube (the same as in 
[Fig. Hfc)] nested in the silver carcass, which is a kind of 
3D wire medium. The metallic structures are embedded 
in silica. 

As the effective refractive index of the split cube in 
carcass does not depend on the slab thickness, it is 
not surprising that the NRW method and our approach 
give results coinciding much better than for the fishnet 
[see Fig. [H^a)]. Nevertheless, effective impedances, and 
therefore permittivities and permeabilities provided by 
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FIG. 8: (Color online). Effective parameters of the split 
cube in carcass negative index MM: refractive index (a), 
impedance (b), permittivity (c) and permeability (d), real 
(black) and imaginary (green/grey) parts. Results by the 
volume-averaged (circles) and surface averaged (triangles) ap- 
proaches are compared with the NRW method (solid lines). 



all three approaches are different [see Figs. [5Jb,c,d)]. 

We should also note that in both cases in the frequency 
ranges beyond the resonances the volume averaging ap- 
proach produces physically sound results. As an illustra- 
tion we note that diamagnetism observed in the 
does not remain in the 3?(/Um)j which is close to conven- 
tional 1 below the resonant region. 



IV. DISCUSSION AND CONCLUSIONS 

We have suggested a novel approach for the extrac- 
tion of effective parameters of metamaterials based on 
the study of dispersion properties of the Bloch waves 
propagating in quasi-periodic structured materials. In 
all the cases with single-mode propagation our approach 
provides solid results for the effective refractive indices, 
which can be attributed to the bulk refractive indices 



of the metamaterials irrespectively of their anisotropy 
and spatial dispersion. Our spectral analysis approach 
is able to retrieve refractive indices for a wide range of 
materials and structure geometries, which can be lossy 
or lossless, dispersive, possess negative permittivity, per- 
meability and refractive index values. The method is 
simple and unambiguous, free from the "branch" and 
Fabry-Perot problems, which are the issues for the re- 
flection/transmission based NRW method. The results 
provided by the NRW method are identical to the results 
obtained by our method in all considered cases except 
for the case of the fishnet MM, where EPs experience 
poor convergence to the bulk values. The single-mode 
propagation of a MM can be checked during the retrieval 
process from the fields mismatch monitoring procedure. 

The spectral analysis serves as a platform for fur- 
ther advance in retrieving EPs. Impedance retrieving 
is very sensitive to the conditions of restoration and can 
lead either to WEP or MEP. Employing surface aver- 
aged fields of the dominating Bloch mode, we obtain 
WEP, which are nearly identical for those retrieved by 
the NRW method, but free from spurious resonances ap- 
pearing from the Fabry-Perot effects in slabs. All what is 
needed for the MEP retrieval accordingly to Ref . pll38j is 
the volume averaging of the electric and magnetic fields 
over the unit cell. Both retrievals (wave and material 
EPs) are performed within a single computational cycle, 
because fields on the unit cells entrance facets or in its 
volumes are available, and they can be exported from 
Maxwell's solver arrays. The approach works for MM 
slabs with thicknesses at least 3-4 monolayers. Our ap- 
proach adequately reveals the typical non-magnetic be- 
havior of metamaterials away from the resonance regions, 
which is problematic for the NRW method. Therefore, 
we anticipate that the proposed approach will become 
a useful tool for the characterization of both wave and 
material effective properties of MMs. 

It should be noted that the magnetic microfields re- 
turned by Maxwells solvers are b//xrj-fields, while the 
volume averaged magnetic field Hva must be used. Pos- 
sible implications of ignoring this fact can be illustrated 
through the Poynting vector calculations, Fig.©. Here 
the fishnet structure from the Subsection 3D is used. 
Poynting vectors are calculated accordingly to three for- 
mulas: 



S,i=9t([ [exh*]dV), 
Jv 

S Z 2 = $t[EvA X Hg^J, 
Sz3 = 5?[EyA X Hy A ]. 



(28) 



Straightforward calculations of the Poynting vector 
give us the negative z-component S Z 3 (see orange line 
with triangles in Fig.©, which means that vectors k and 
S are parallel in the negative index domain. This is ex- 
actly what can happen if the wrong formulation of the 
Poynting vector through vector b is used as it is pointed 
in Refs . |28lJ32j | . The consequences of this is not only the 
wrong direction of the flux, but also the negative value 
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FIG. 9: (Color online) z-component of the Poynting vector 
of the fishnet negative-index MM: volume averaged Poynting 
vector (red line with circles), correctly defined Poynting vector 
for the fundamental Bloch harmonic (black line with squares) 
and flux calculated through the volume averaged electric and 
magnetic fields of the fundamental Bloch harmonic (orange 
line with triangles). 



of the 5i(z), because flux and impedance are connected 
through the expression 



S z3 = Vl(e z [EvA x H^]) = 

= ft(E VA H^ A ) = Z n(z w )\H VA \ 



(29) 



However, employment of the volume averaged electric 
and surface averaged magnetic fields improves the situa- 
tion (black line with squares). The Poynting vector S Z 2 
calculated through them is very close to the averaged 
microscopic flux S zl (red line with circles). Such calcula- 
tions confirm the fact that on the grid level, microfields 
b and h differ only by a constant. But fields averaged 
over a macrovolume bear principle differences. 

The most intriguing part is the direct comparison be- 
tween effective parameters restored with formulas (I12I) . 
(fl~3f and ([22]). In Fig. ([10]) we plot results for three dif- 
ferent cases of impedance restoration and include also 
the NRW restoration data. In fact, the volume-averaged 
fields provide the incorrect result (stars), with nega- 
tive 3?(z) and double anti-resonances in 9(e) and Q(fJ.). 
The situation is improved when the surface-averaged 
(transverse-averaged) fields are taken (triangles) instead 
of bulk fields in concordance with the finding in™. There 
is still one faint "attempt" of an antiresonance with de- 
creasing of 3(e). And there is completely no antireso- 
nance, when using the formula (I21[) . The corresponding 
curves are designated by circles in Fig.tfTT)]). 

Unfortunately this approach cannot be accepted as 
universal retrieving method, because in some cases (see 
Split Cube case in Section 3C) it fails. More deep anal- 
ysis in the failure of formula (|2"2"|) is needed, but it lies 
beyond the scope of this paper. 
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FIG. 10: (Color online) Effective parameters of the split 
cube in carcass negative index MM: impedance (a), permit- 
tivity (b) and permeability (c), real (black) and imaginary 
(green/grey) parts. Results are obtained by formula (|22l) (cir- 
cles), JT2) (triangles), ([13]) (stars) and the NRW method (solid 
lines)approaches . 



We should admit that a direct extension of our ap- 
proach for the experimental characterization of MMs in 
the optical range is challenging, since there are no such 
small electric and magnetic field detectors that could be 
placed inside the MM unit cell without noticeable influ- 
ence on its functionality. Nevertheless, as the radio and 
microwave frequency range, it is possible to record the 
fields at the spatial points inside the metamaterial^i, en- 
abling the direct application of our approach. 
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